subroutine readinput
    use allocs
    use define
    implicit none

    include "pspec.FI"
    include "files.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "hmgz.FI"
    include "mcbase.FI"
    include "stsol.FI"
    include "CMFD_LARSEN.FI"
    
    integer unit, i, j, k, index, index_cmfd, g, ichar1st, temp, ratio, base, left, right
    character*200 logstring
    character*2048 oneline, cardname
    character(1) probe, sc(2048)
    character(1) DOT,BANG,BLANK,SLASH,AST,POUND
    parameter (DOT='.',BANG='!',BLANK=' ',SLASH='/',AST='*')
    parameter (POUND='#')
    integer :: pin(200), assembly(200,200), id, val, npasm(200)
    integer :: core(1000), nasm, ndataf
    integer :: totasm, nmpin, totpin
    real(8) :: pinsize
    integer, parameter :: fref = 90
    integer :: i_c
    real(8) :: cphi
    equivalence (probe,oneline)
    equivalence (sc(1),oneline)

    real(8), pointer, dimension(:,:) :: tsigth
    real(8), pointer, dimension(:,:,:) ::tt
    integer :: ci
    
    localfn = ' '
    
!    call getarg(1, localfn)
!    ichar1st=ichar(localfn(1:1))
!    if(ichar1st.eq.0 .or. ichar1st.eq.32) localfn='input.txt'
    localfn='input.txt'
    
    unit =FIinp    
    unifsrc = .false.
    open(unit, file=localfn,status='old')
    read(unit,*) FDBIACT, FDBACT, PDFREF, WRTREF, delkw, FLTIGS, SUPERHST
    read(unit,*) SRCSPLT, FMSHAPE, NMVWND, STMVWND, cmfdtype, delay_cmfd
    read(unit,*) FXDSRC, RSTDHAT, MOREPTC, NMPTCR, NMCYCLE
    read(unit,*) pinsize, nmpin
    read(unit,*) albedo  (1:2)
    read(unit,*) cck, ccpsi
    read(unit,*) history, inactive, active
    if (cmfdtype>0) lgcmfd = .true.
    
    !if (SRCSPLT <1 .or. SRCSPLT>2) SRCSPLT = 1
    if (SRCSPLT <1) SRCSPLT = 1
    npasm = 0
    g_splt = .false.
    gcmfd_opts = 0
    impl_cpt = .true.
    do while (TRUE)
        read(unit,'(a2048)') oneline
        if(probe.eq.DOT .or. probe.eq.SLASH) exit
        if(probe.eq.BANG .or. oneline.eq.BLANK .or. probe.eq.POUND) cycle
        
        read(oneline, *) cardname
        call toupper1(cardname)
        select case(cardname)
            case('PIN')
                read(oneline, *) cardname, id, val
                pin(id) = val
            case('ASSEMBLY')
                ndataf=nfields(oneline)-2
                read(oneline, *) cardname, id, assembly(1:ndataf, id)  
                npasm(id) = ndataf
!                if (npasm .ne. 0 .and. npasm .ne. ndataf) write(*,*) "Geometry error : assembly ", id
!                npasm = ndataf
            case('CORE')
                totasm=nfields(oneline)-1
                read(oneline, *) cardname, core(1:totasm)
            case('UNIFORMFSRC')
                unifsrc = .true.
            case('GSPLIT')
                g_splt = .true.
            case('GCMFD')
                ndataf=nfields(oneline)-1
                read(oneline, *) cardname, gcmfd_opts(1:ndataf)
            case('IMPLICIT_CAPTURE')
                read(oneline, *) cardname, impl_cpt
        end select
    enddo
    close(unit)
 
    
    call dmalloc(hx,      totasm)
    call dmalloc(dh,      totasm)
    call dmalloc(nmesh,   totasm)
    call dmalloc(comp,    totasm)
    call dmalloc(nm_cmfd, totasm)
    call dmalloc(dh_cmfd, totasm)
    
    
    totpin = 0
    do i=1, totasm
        if (npasm(core(i)) .eq. 0) write(*,*) "Error in assembly input sheet! "
        totpin = totpin+npasm(core(i))
    end do
    totmesh = totpin*nmpin
    totm_cmfd = totasm


    ! flat fission 
    allocate(pwadj(totmesh), pwadjr(totmesh), nfn(totmesh))
!        pwadj(1:10)     =1._8
!        pwadj(11:20)    =1._8
!        pwadj(21:120)  =1._8
!        pwadj(121:130) =1._8
!        pwadj(131:140) =1._8
        pwadj=1._8
!    pwadj=pwadj*totmesh/sum(pwadj)
        popctrl=1._8*totmesh/sum(pwadj)
    do i=1, totmesh
        pwadjr(i)=1._8/pwadj(i)
    enddo
    ! normalize
    pwadjr=pwadjr*totmesh/sum(pwadjr)
    do i=1, totmesh
        pwadj(i)=1._8/pwadjr(i)
    enddo
    !pwadj=pwadj*sum(pwadj)/totmesh

    
    rcf = totmesh/totm_cmfd
    
    write(*,*) "# of fine mesh   : ", totmesh
    write(*,*) "# of coarse mesh : ", totm_cmfd
    
    
    sizeprb = pinsize*totpin
    dh      = pinsize/nmpin
    do i=1, totasm
        dh_cmfd(i) = pinsize*npasm(core(i))
    enddo
    
    call dmalloc0(complist, 0,totmesh+1)
    call dmalloc (h, totmesh)
    call dmalloc0(hac, 0, totmesh)
    call dmalloc (hfmcmid, totmesh)
    call dmalloc (g_av, totmesh+1)
    call dmalloc(hfmcorg, totmesh)

!    call dmalloc (compl_cmfd, totm_cmfd)
    call dmalloc (h_cmfd, totm_cmfd)
    call dmalloc0(hac_cmfd, 0, totm_cmfd)

    
    call dmalloc(psi_f, totmesh)
    call dmalloc(psi_f2, totmesh)
    call dmalloc(psifAC, totmesh)
    call dmalloc(psi_c, totm_cmfd)
    call dmalloc(psi_cp, totm_cmfd)
   
    index = 1
    hac(0) = 0.
    do i=1, totasm
        do j=1, npasm(core(i))
            do k=1, nmpin 
                complist(index) =pin(assembly(j, core(i)))
                h(index) = dh(i)
                hac(index) =hac(index-1) + h(index)
                index = index +1
            end do
        end do 
    end do
!    hfmcmid = h/2



    
    hac_cmfd(0)= 0.
    do i=1, totasm
        h_cmfd(i) = dh_cmfd(i)
        hac_cmfd(i) = hac_cmfd(i-1) + h_cmfd(i)
    end do
    
    g_av(1) = 2._8/h_cmfd(1)
    do i=2, totm_cmfd
        g_av(i) = 2._8/(h_cmfd(i-1)+h_cmfd(i))
    end do
    g_av(i) = 2._8/h_cmfd(i-1)
    
!    g_av(1) = 2/h(1)
!    do i=2, totmesh
!        g_av(i) = 2/(h(i-1)+h(i))
!    end do
!    g_av(i) = 2/h(i-1)


    do i=1, totm_cmfd
        hfmcorg(rcf*(i-1)+1:i*rcf) = hac_cmfd(i-1)
        hfmcmid(rcf*(i-1)+1:i*rcf) = h_cmfd(i)/2
    enddo
        
    unit =FXsec
    open(unit, file='xsec.txt',status='old')
    read(unit,*) ng, nComp
    
    call dmalloc(gtally_arr, 3, inactive+active)


    call dmalloc0(Xdif, 0,nComp+1, 1,ng)
    call dmalloc0(Xrmv, 0,nComp+1, 1,ng)
    call dmalloc0(Xnuf, 0,nComp+1, 1,ng)
    call dmalloc0(Xchi, 0,nComp+1, 1,ng)
    call dmalloc0(Xsct, 0,nComp+1, 1,ng)
    call dmalloc0(Xtot, 0,nComp+1, 1,ng)
    call dmalloc0(Xabs, 0,nComp+1, 1,ng)
    call dmalloc0(XsctM, 1,ng, 1,ng, 0,nComp+1)
    call dmalloc0(keffg, 0,nComp+1, 1,ng)
    
    do i=1 , nComp
        read(unit,*)
        do g=1, ng
            read(unit,*) temp, Xdif(i,g), Xrmv(i,g), Xnuf(i,g), Xchi(i,g)
            keffg(i,g) = Xnuf(i,g)/3/Xdif(i,g)
        end do
        do g=1, ng
            read(unit,*) temp, XsctM(1:ng, g, i)
            XsctM(g,g,i) = 1./3/Xdif(i,g)-Xrmv(i,g)
            Xsct (i, g) = sum(XsctM(1:ng, g, i))
        end do
    end do
    close(unit)
    Xtot(1:ncomp,:) = 1./(3*Xdif(1:ncomp,:))
    Xabs = Xtot-Xsct
    
    call dmalloc(mt_f, 3,   totmesh,   ng)
    call dmalloc(jf,        totmesh*2, ng)
    call dmalloc(pjf,       totmesh*2, ng)
    call dmalloc(Dtilde_f,  totmesh*2, ng)
    call dmalloc(phi_f,     totmesh,   ng)
    call dmalloc(phi_fnm,   totmesh,   ng)
    call dmalloc(phi_f2,    totmesh,   ng)

    call dmalloc(adjp,      totmesh)
    call dmalloc(srcPDF,    totm_cmfd)

    call dmalloc(trk_est,   totmesh,   ng)
#ifdef _2NDHRMN
    call dmalloc(trk_est2,  totmesh,   ng)
    call dmalloc(phi_f2nd,  totmesh,   ng)
#endif
!    call dmalloc(trk_estT,  totmesh,   ng)
    
    call dmalloc(mt_cmfd,3, totm_cmfd,   ng)
    call dmalloc(Dtilde_c,  totm_cmfd*2, ng)
    call dmalloc(Dhat,      totm_cmfd*2, ng)

    call dmalloc(mtfull, totm_cmfd*ng, totm_cmfd*ng)
    call dmalloc(mtflu, totm_cmfd*ng, totm_cmfd*ng)

! Variables for FMC
    call dmalloc(phi_s,  totmesh+1, ng)
    call dmalloc(phi_sa, totm_cmfd+1, ng)
    call dmalloc(mmt1st, totmesh+1, ng)
    call dmalloc(mmt2nd, totmesh+1, ng)
    call dmalloc(mmtAbsr, totmesh+1, ng)
    call dmalloc(mmtNuF,  totmesh+1, ng)
    call dmalloc(adj_fmc, totmesh+1)

!! Variables for DF
!    call dmalloc(coef_df, totm_cmfd*2, ng)

    cntacc = 0
    
    allocate(ntrn(history*3, 3));
    
    do g=1, ng
        do i=1, totmesh
        phi_f(i,g) = 1.
        end do
    end do
    
    
    call dmalloc(jcvt2cmfd, totmesh*2)
    jcvt2cmfd = FALSE
    left = 1
    do i=1, totasm
!        npasm(core(i))
!        ratio = npasm*nmpin*2
!        left  = (i-1)*ratio+1
!        right = i*ratio
        right = left + nmpin*npasm(core(i))*2-1
        jcvt2cmfd(left ) = TRUE
        jcvt2cmfd(right) = TRUE
        left = right+1
    end do
    
    call dmalloc0(fi2co, 0,totmesh+1)
    index =1
    index_cmfd = 1
!    ratio = npasm*nmpin
    do i=1, totasm
        fi2co(index: index+ nmpin*npasm(core(i))-1) = index_cmfd
        index = index+ nmpin*npasm(core(i))
        index_cmfd = index_cmfd+1
    end do
     
    ! for test MC
    call dmalloc(absn, totmesh, ng)
    call dmalloc(srcn, totmesh, ng)
    call dmalloc(sctn, ng, ng, totmesh)
    


    if (PDFREF) then
        logstring = "Read source distribution from reference solution"
        call log_(logstring)    
        open(fref, file='REF.txt', status='old')
        call dmalloc(refphi, totmesh, ng)
        call dmalloc(refpsi, totmesh)
        do g = 1, ng
            read(fref, *) refphi(:,g)
        end do
        close(fref)

        call dmalloc(refc, totm_cmfd, ng)
        call dmalloc(refPDF, totm_cmfd)
        refc = 0._8
        refPDF = 0._8
        refpsi = 0._8
        do g=1, ng
            do i=1, totmesh
                refpsi(i) = refpsi(i) + Xnuf(complist(i), g) * refphi(i, g) * h(i)
                i_c = fi2co(i)
                cphi  = refphi(i,g)*h(i)/h_cmfd(i_c)
                ! cycle flux are used to obtain initial guess
                refc(i_c, g) = refc(i_c, g) + cphi
                refPDF(i_c) = refPDF(i_c) + xnuf(complist(i),g)*cphi
            end do
            refpsi = refpsi/sum(refpsi)
            refPDF = refPDF/sum(refPDF)
!            do i=1, totm_cmfd
!                refPDF(i) = refPDF(i-1)+refPDF(i)
!            end do
        end do

        l2ref_f = 0
        do g=1, ng
            do i=1, totm_cmfd
                l2ref_f = l2ref_f + refphi(i, g)**2
            end do
        end do
        l2ref_f = sqrt(l2ref_f)
        sumref_f = sum(refphi)
    
    endif   
    
    call dmalloc(refflux, totmesh, ng)
!    open(fref, file='REFMC.txt', status='old')
    do g=1, ng
!        read(fref, *) refflux(:,g)
    enddo
    
    refflux=refflux/sum(refflux)
    ! homogenization
    allocate(hpsi(totm_cmfd))
    allocate(hphi(totm_cmfd, ng))
    allocate(hphi2(totm_cmfd, ng))
    allocate(hphis(totm_cmfd*2, ng))
    allocate(hj(totm_cmfd*2, ng))
    allocate(hpj(0:totm_cmfd*2+1, ng))
    allocate(hDif(totm_cmfd, ng))
    allocate(hTot(totm_cmfd, ng))
    allocate(hRmv(totm_cmfd, ng))
    allocate(hNuF(totm_cmfd, ng))
    allocate(hmSct(ng,ng, totm_cmfd))
    allocate(hTot_volavg(totm_cmfd, ng))
    
    ! accumulation
    allocate(psifAC(totmesh))
    allocate(phifAC(totmesh, ng))
    allocate(phif2AC(totmesh, ng))
    allocate(jfAC(totmesh*2, ng))
    allocate(pjAC(totmesh*2, ng))
    allocate(phisAC(totm_cmfd+1, ng))
    allocate(phisaAC(totm_cmfd+1, ng))
    allocate(mmt1stAC(totmesh+1, ng))
    allocate(mmt2ndAC(totmesh+1, ng))
    allocate(mmtAbsrAC(totmesh+1, ng))
    allocate(mmtNuFAC(totmesh+1, ng))
    
    if (RSTDHAT) then
        allocate(phird(totmesh,ng,inactive))
        allocate(jrd(totmesh*2,ng,inactive))
    endif
        
    call dmalloc0(stmc, 0,totmesh, 1,ng)
    call dmalloc0(stsqmc, 0,totmesh, 1,ng)

    call dmalloc0(stmcf, 0,totmesh)
    call dmalloc0(stsqmcf, 0,totmesh)

    call dmalloc0(stcmfd, 0,totm_cmfd, 1,ng)
    call dmalloc0(stsqcmfd, 0,totm_cmfd, 1,ng)
    
    call dmalloc0(stcmfdf, 0,totm_cmfd)
    call dmalloc0(stsqcmfdf, 0,totm_cmfd)

    call dmalloc0(stfmc, 0,totmesh+1)
    call dmalloc0(stsqfmc, 0,totmesh+1)
    
    
#ifdef CMFDLARSEN
    call dmalloc0(dtildeL, 0, totm_cmfd, 1,ng)
    call dmalloc0(dhatL,   0, totm_cmfd, 1,ng)
    call dmalloc0(phi2Ls,  0, totmesh,   1,ng)
    call dmalloc0(phi1Ls,  0, totmesh,   1,ng)
    call dmalloc0(phi2Lsh, 0, totm_cmfd, 1,ng)
    call dmalloc(phi1La,      totmesh,     ng)
    call dmalloc(phi1tLa,     totmesh,     ng)
    call dmalloc(phi1Lah,     totm_cmfd,   ng)
    call dmalloc(phi1Laacc,     totm_cmfd,   ng)
    call dmalloc0(delx_2,  0, totmesh)
    call dmalloc(htjr,     totmesh,     ng)
    
    
    delx_2(0)=h_cmfd(1)/2
    delx_2(totm_cmfd)=h_cmfd(totm_cmfd)/2
    do i=1, totm_cmfd-1
        delx_2(i) = (h_cmfd(i)+h_cmfd(i+1))/2
    enddo
#endif    
!   GCMFD with tent function
    call dmalloc0(sigth_hf, 0,totm_cmfd, 1,ng)
    call dmalloc0(fphi1_hf, 1,2, 0,totmesh, 1,ng)
    call dmalloc (gphi0_hf, 3, totmesh, ng)
    call dmalloc (phi0, totm_cmfd, ng)
    call dmalloc (phi2, totm_cmfd, ng)
    call dmalloc (phi2_hf, 2,totm_cmfd, ng)
    call dmalloc (phi1p, 2, ng)
    call dmalloc (hcmfdhf, totm_cmfd)
    call dmalloc0(phi1c, 0,totm_cmfd, 1,ng)
    
    call dmalloc0(fphi2fs, 1,2, 0,totmesh, 1,ng)
    call dmalloc(phi2f, totmesh, ng)
    
    do i=1, totm_cmfd
        hcmfdhf(i)= hac_cmfd(i)-h_cmfd(i)/2
    enddo
    
    
    call dmalloc0(h_hf,    0,totm_cmfd)
    call dmalloc0(tt, 1,2, 0,totm_cmfd, 1,ng)
!    gphi0_hf(1,1,1)=10._8

    h_hf(0)=h_cmfd(1)/2
    do i=1, totm_cmfd-1
        h_hf(i)=(h_cmfd(i)+h_cmfd(i+1))/2
    enddo
    h_hf(totm_cmfd)=h_cmfd(totm_cmfd)/2
    
    allocate(tsigth(totm_cmfd,ng))
    tsigth = 0
    do g=1, ng
    do i=1, totmesh
        ci = fi2co(i)
        tsigth(ci,g)=tsigth(ci,g)+xtot(complist(i), g)*h(i)
    enddo
    enddo

    do g=1, ng
        sigth_hf(0,g)=(tsigth(1,g))/2        
        do i=1, totm_cmfd-1
            sigth_hf(i,g)=(tsigth(i,g)+tsigth(i+1,g))/2
        enddo
        sigth_hf(totm_cmfd,g)=(tsigth(totm_cmfd,g))/2        
    enddo
    deallocate(tsigth)
    
    call dmalloc(xsum,totm_cmfd, ng)
    call dmalloc(xsqsum,totm_cmfd, ng)
    call dmalloc(xysum,totm_cmfd*2, ng)
    call dmalloc(ysum,totm_cmfd*2, ng)

    call dmalloc(beta,totm_cmfd*2, ng)

contains

    subroutine toupper1(aa)
    !
    ! convert lowercase string to upZpercase
    integer:: INDXA, IDNXZ, lenaa, ia
    parameter (INDXA=97,IDNXZ=122)
    character aa*(*)
    !
    lenaa=len_trim(aa)
    i=1
    do while (aa(i:i).ne.' ' .and. i.le.lenaa)
     ia=ichar(aa(i:i))
     if(ia.ge.INDXA) aa(i:i)=char(ia-32)
     i=i+1
     if(i.gt.lenaa) return
    enddo
    end subroutine

    integer function nfields(aline)
    character(2048)  aline
    logical nonblankd,nonblank,multidata
    integer ncol, n, multidcol, nmult
    !
    nonblankd=.false.
    multidata=.false.
    oneline=aline
    ncol=len_trim(aline)
    n=0
    do i=ncol,1,-1
     if(sc(i).eq.BLANK .or. ichar(sc(i)).eq.9) then !ichar(tab)=9
        nonblank=.false.
     else
        if(sc(i).eq.AST) then
           multidata=.true.
           multidcol=i
        endif
        nonblank=.true.
        if(sc(i).eq.BANG) then
           n=-1
           nonblankd=.true.
        endif
     endif
     if((.not.nonblankd.and.nonblank) .or. &
          (nonblankd.and..not.nonblank)) then
        n=n+1
        if(multidata.and.(nonblankd.and. .not.nonblank)) then
           read(oneline(i+1:multidcol-1),*) nmult
           n=n+(nmult-1)*2
           multidata=.false.
        endif
     endif   
     nonblankd=nonblank
    enddo
    if(mod(n,2).ne.0) then
    nfields=n/2+1
    else
    nfields=n/2
    endif
    !
    return
    end function
!logical function ifnumeric(aline)
!character aline*(*),al*256
!integer ::iascii
!!
!ifnumeric=.false.
!al=aline
!do i=1,100
! iascii=ichar(al(i:i))
! if(al(i:i).ne.BLANK .and. iascii.ne.9 .and. iascii.ne.45 ) then  !determine if the first character is numeric
!    if((iascii-48)*(iascii-57) .le. 0) ifnumeric=.true.
!    return
! endif
!enddo
!!
!return
!end function
  
end subroutine